Os índigenas Pima - diabetes

Gabriel de Jesus Pereira

2025-02-25

Sobre o banco de dados

Sobre o banco de dados

  • Foi originado do National Institute of Diabetes and Digestive and Kidney Diseases.

  • Essa base foi criada com base nos Pima, um grupo de nativos americanos que vivem em uma área que atualmente abrange o centro e o sul do estado do Arizona.

  • A base possui 768 observações e 9 colunas. O objetivo é classificar se um paciente tem diabetes com base em algumas características do paciente.

Descrição das colunas

  • Pregnancies: Número de gestações da paciente.

  • Glucose: Concentração de glicose plasmática em teste oral de tolerância à glicose.

  • BloodPressure: Pressão arterial diastólica (mm Hg).

  • SkinThickness: Espessura da dobra cutânea do tríceps (mm).

  • Insulin: Nível sérico de insulina (mu U/ml).

  • BMI: Índice de Massa Corporal (peso em kg / altura em m², IMC).

  • DiabetesPedigreeFunction: Histórico familiar de diabetes (mede a predisposição genética).

  • Age: Idade do paciente.

  • Outcome: 0 é não diabético e 1 é diabético.

Análise exploratória de dados

Quantidade de valores ausentes

library(dplyr)
library(ggplot2)
library(kknn)
library(naivebayes)
library(glmnet)
library(discrim)
library(tidyverse)
library(purrr)
library(tidymodels)
library(visdat)
library(corrplot)

data <- read_csv("diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
cols <- c("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")
data[cols][data[cols] == 0] <- NA


Ausentes <- data |>
  vis_dat() +
  labs(
    y = "Observações"
    ) +
  scale_fill_manual(
    values = viridis::mako(10)[5:6],
    labels = c("Qualitativa", "Quantitativa", "Ausentes")
    ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
  )

Quantidade de valores ausentes

Apenas as colunas de glicose, pressão arterial, espessura da dobra cutânea do tríceps, nível de insulina e IMC. Aquela que mais possui valores ausentes é a variável de nível de insulina.

Frequência

set.seed(123)
data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) %>%
  step_impute_median(all_numeric_predictors())

prep_rec <- prep(rec)
train_data_eda <- bake(prep_rec, new_data = train_data)
test_data_eda <- bake(prep_rec, new_data = test_data)

Frequência <- train_data_eda |>
  ggplot(aes(x = forcats::fct_infreq(as.factor(Outcome)))) +
  geom_bar(
        aes(
          y = after_stat(count) / sum(after_stat(count)),
          fill = Outcome
          ),
        show.legend = FALSE,
  ) +
  theme_bw() +
  theme(axis.title.x = element_blank()) +
  scale_fill_manual(
    values = viridis::mako(5)[3:5]
    ) +
  labs(y = NULL, x = NULL)

Frequência

A base de dados possui um problema de desbalanceamento das classes. A classe com maior frequência é de não diabéticos, com pouco mais de 60%.

Distribuição

Distribuição <- train_data_eda |>
  pivot_longer(
    !Outcome,
    names_to = "variable",
    values_to = "values"
  ) |>
  mutate(Outcome = as.factor(Outcome)) |>
  ggplot(
    aes(
      x = values,
      fill = Outcome
      )
    ) +
  geom_violin(
    aes(y = variable),
    draw_quantiles = c(0.25, 0.5, 0.75),
    ) +
  labs(x = "", y = "") +
  facet_wrap(vars(variable), scales = "free") +
  scale_fill_manual(
    values = viridis::mako(5)[3:5],
    name = "Outcome"
    ) +
  theme_bw() +
  theme(axis.text.y = element_blank())

Distribuição

As variáveis apresentam uma assimetria negativa, sendo a variável de pressão arterial a única mais simétrica. Além disso, observando a distribuição dos diabéticos e não diabéticos para cada variável, eles apresentam um padrão de comportamento bastante similar.

Correlação

data_cor <- train_data_eda |>
  (\(x) {
     scaled <- x |>
       select_if(is.numeric) |>
       scale() |>
       as_tibble()

     x |>
       select_if(is.character) |>
       cbind(scaled)
  })()

data_cor |>
  select_if(is.numeric) |>
  cor() |>
  corrplot(
      method="circle",
      col = viridis::mako(200),
      diag = FALSE,
      addgrid.col = "#00000040",
      tl.pos = "l",
      tl.cex = 0.6,
      tl.col = "black",
      tl.offset = 0.25,
      cl.length = 9,
      cl.cex = 0.6,
      cl.ratio = 0.25,
      family = "mono"
  )

Correlação

A maioria das variáveis independentes não apresentam uma alta correlação. No entanto, os níveis de insulina tem uma alta correlação com a glicose, o IMC tem uma alta correlação com a espessura da dobra cutânea do tríceps e a idade tem uma alta correlação com a quantidade de gravidez.

Métricas utilizadas para avalização dos modelos

Métricas utilizadas para avalização dos modelos

  • ROC AUC

  • Acurácia

  • Brier Score: \(BS = \frac{1}{N}\sum_{t=1}^N\left(f_t - o_t\right)^2\). \(N\) é o tamanho da amostra classificada, \(f_t\) é a probabilidade predita e \(o_t\) é o valor observado.

Preparação dos dados

Preparação dos dados

data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) |>
  step_impute_median(all_numeric_predictors()) |>
  step_mutate(across(all_numeric_predictors(), log1p)) |>
  step_mutate(
    N0 = BMI * SkinThickness,
    N8 = Age / Insulin,
    N13 = Glucose / DiabetesPedigreeFunction,
  ) |>
  step_nzv(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())

Preparação dos dados

  • Os dados foram divididos entre teste e treino, com 20% para o conjunto de teste e estratificado pela variável dependente.

  • Os valores ausentes foram substituídos pela mediana das variáveis.

  • Foram criadas novas variáveis, idade por nível de insulina, glicose por predisposição de ter diabetes e o produto entre IMC e a espessura da dobra cutânea do tríceps.

  • Todas as variáveis foram transformadas com \(\log \left(1 + x\right)\) e normalizadas. Além disso, foi utilizado o utilizado o step_nzv, que serve para remover variáveis com variância próximo de zero.

Algoritmos e validação cruzada

cv_folds <- vfold_cv(train_data,
                     v = 10,
                     strata = Outcome
                     )

knn_spec <- nearest_neighbor(neighbors = tune(), weight_func = tune()) |>
  set_engine("kknn") |>
  set_mode("classification")

nbayes_spec <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
  set_engine("naivebayes") |>
  set_mode("classification")

lda_spec <- discrim_linear() |>
  set_engine("MASS") |>
  set_mode("classification")

qda_spec <- discrim_quad() |>
  set_engine("MASS") |>
  set_mode("classification")

logistic_spec <- logistic_reg(
  engine = "glmnet",
  penalty = tune(),
  mixture = tune()
)

wf = workflow_set(
  preproc = list(rec),
  models = list(
    knn_fit = knn_spec,
    nbayes_fit = nbayes_spec,
    lda_fit = lda_spec,
    qda_fit = qda_spec,
    linear_fit = logistic_spec
    )
  )  |>
  mutate(wflow_id = gsub("(recipe_)", "", wflow_id))

grid_ctrl = control_grid(
  save_pred = TRUE,
  parallel_over = "resamples",
  save_workflow = TRUE
)

grid_results = wf  |>
  workflow_map(
    seed = 123,
    resamples = cv_folds,
    grid = 50,
    control = grid_ctrl
  )
autoplot(grid_results, metric = "roc_auc") +
  theme_bw() +
  labs(y = "Métrica", x = "")
autoplot(grid_results, select_best = TRUE) +
  theme_bw() +
  labs(y = "Métrica", x = "")

Algoritmos e validação cruzada

  • A validação cruzada foi realizada com 10 folds, estratificado pela variável dependente e foi criado um grid de 50 combinações entre os hiperparâmetros. Assim, foram otimizados os algoritmos a quantidade de vizinhos e o tipo de distância utilizada. Foram otimizados, para o naive bayes, o parâmetro de suavidade para o limite da classe e a correção da suavidade chamado de Laplace. Foram também otimizados os hiperparâmetros de regularização da regressão logística.

Resultado geral da otimização

Resultado na base de teste

results_acc = workflowsets::rank_results(
  grid_results,
  select_best = TRUE,
  rank_metric = "roc_auc"
  ) |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, mean, std_err, model, rank)

model_ids <- c("linear_fit", "knn_fit", "nbayes_fit", "lda_fit", "qda_fit")

best_sets <- map(model_ids, ~ grid_results %>%
                   extract_workflow_set_result(.x) %>%
                   select_best(metric = "roc_auc"))

names(best_sets) <- model_ids

test_results <- function(rc_rslts, fit_obj, par_set, split_obj) {
  rc_rslts %>%
    extract_workflow(fit_obj) %>%
    finalize_workflow(par_set) %>%
    last_fit(
      split = split_obj,
      metrics = metric_set(accuracy, roc_auc, brier_class
      )
    )
}

test_results_list <- map2(model_ids, best_sets,
                          ~ test_results(grid_results, .x, .y, data_split))

metrics_table <- map_dfr(test_results_list, collect_metrics, .id = "model_id") |>
  dplyr::select(model_id, .metric, .estimate) |>
  pivot_wider(names_from = .metric, values_from = .estimate) |>
  mutate(across(where(is.numeric), round, 4)) |>
  mutate(Modelo = c("Logistic", "K-NN", "Naive Bayes", "LDA", "QDA")) |>
  relocate(Modelo) |>
  dplyr::select(-model_id)

colnames(metrics_table) <- c("Modelo", "Accuracy", "ROC AUC", "Brier Score")

metrics_table <- as_tibble(metrics_table)
metrics_table |>
  arrange(desc(`ROC AUC`)) |>
  knitr::kable()
Modelo Accuracy ROC AUC Brier Score
LDA 0.7857 0.8652 0.1443
Logistic 0.7922 0.8598 0.1483
K-NN 0.7987 0.8451 0.1545
QDA 0.7727 0.8238 0.1663
Naive Bayes 0.7013 0.8128 0.1842

Resultado na base de teste

O modelo que obteve o melhor resultado foi o discriminante linear, com base na curva ROC e o brier. Na base de teste, ele obteve um ROC AUC de 0,8652, a métrica de brier 0,1443 e obteve uma acurácia de 78,57%. O que pior perfomou foi o Naive Bayes, obtendo as piores estatísticas para todas as métricas consideradas.

Hiperparâmetros selecionados dos modelos tunados

# Hiperparâmetros para o K-NN

grid_results |>
  extract_workflow_set_result("knn_fit") |>
  select_by_one_std_err(desc(neighbors), metric = "roc_auc") |>
  knitr::kable()
neighbors weight_func .config
15 rank Preprocessor1_Model30
# Hiperparâmetros do naive bayes

grid_results |>
  extract_workflow_set_result("nbayes_fit") |>
  select_by_one_std_err(desc(Laplace), metric = "roc_auc") |>
  knitr::kable()
smoothness Laplace .config
0.9762805 2.978423 Preprocessor1_Model23
# Hiperparâmetros da regressão logística com penalização

grid_results |>
  extract_workflow_set_result("linear_fit") |>
  select_by_one_std_err(desc(mixture), metric = "roc_auc") |>
  knitr::kable()
penalty mixture .config
5.8e-06 0.9931674 Preprocessor1_Model50

Referência

Pima Indians Diabetes. Baixado através do Kaggle em: https://www.kaggle.com/datasets/uciml/pima-indians-diabetes-database.